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Abstract - We use different numerical approaches to calculate the double occupancy and mag- 
netic susceptibility as a function of a bias voltage in an Anderson impurity model. Specifically, 
we compare results from the Matsubara- voltage quantum Monte-Carlo approach (MV-QMC), the 
scattering-states numerical renormalization group (SNRG), and real-time quantum Monte-Carlo 
(RT-QMC), covering Coulomb repulsions ranging from the weak-coupling well into the strong- 
coupling regime. We observe a distinctly different behavior of the double occupancy and the 
magnetic response. The former measures charge fluctuations and thus only indirectly exhibits the 
Kondo scale, while the latter exhibits structures on the scale of the equilibrium Kondo tempera- 
ture. The Matsubara- voltage approach and the scattering-states numerical renormalization group 
yield consistent values for the magnetic susceptibility in the Kondo limit. On the other hand, 
all three numerical methods produce different results for the behavior of charge fluctuations in 
strongly interacting dots out of equilibrium. 
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Introduction. — The advances in nanostructuring 
of heterogeneous semiconductors and in the handling of 
molecules have made it possible to reproducibly build care- 
fully designed nanometer-scale devices. These generically 
consist of a few locally interacting degrees of freedom, 
for example in a quantum dot in contact with macro- 
scopic leads. The spatial confinement of the quantum 
dot electrons to a few nanometers implies a small elec- 
trical capacitance C and, hence, a sizable charging energy 
U — e^/C. The attraction of such devices stems from 
their highly controllable properties fT|2l . Substantially in- 
creasing the coupling to the leads but still maintaining the 
charge fluctuation scale below charging energy tunes the 
quantum dots into the experimentally accessible Kondo- 
regime (3||6|. This regime is characterized by the lifting 
of the Coulomb blockade [l] in the quantum transport at 
temperatures below a dynamically generated small energy 
scale, called Kondo temperature T^, which ubiquitously 
shows up in physical properties [7j. 
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The experiments performed on mesoscopic systems typ- 
ically are measurements of transport properties in the 
presence of external fields and voltage bias, making a 
theoretical description in terms of non-equilibrium sta- 
tistical physics mandatory. The theoretical challenge of 
the Kondo regime is related to the change of ground 
state [t] upon cooling the system from an intermediate- 
temperature local-moment regime to the low-temperature 
regime which manifests itself in the lifting of the Coulomb 
blockade at zero bias. This crossover cannot be reliably 
accessed by any finite order perturbation theory in the 
Coulomb repulsion and requires more sophisticated an- 
alytical methods such as Bethe ansatz [8 or numerical 
methods such as Wilsons numerical renormalization group 
(NRG) :9 . 

However, no exact solution for the transport proper- 
ties through quantum dots at finite bias exists for models 
of interacting quantum many-body systems out of equi- 
librium. Over the past two decades, several approaches 
have been developed to approximately or numerically solve 
such models, ranging from perturbation theory 10-12 
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and renormalization approaches 13 - 15 to various numcr 



ical techniques 16 -22 . Most of these approaches work 
well in certain limits, where Kondo physics is either not 
yet relevant or ~ due to strong external fields - already 
suppressed. Accessing the crossover regime where exter- 
nal fields, in particular the bias voltage across the dot, are 
of the order of the Kondo temperature, remains a great 
challenge. 

In this paper we present results for static quantities of 
a quantum dot in steady-state non equilibrium for dot 
parameters, temperatures, and voltages that fall precisely 
into this challenging regime. One approach employed here 

It 



23 



has been recently proposed by Han and Heary 
maps the steady-state non-equilibrium system onto an in- 
finite set of auxiliary equilibrium statistical-physics prob- 
lems. The latter are solved by a continuous-time quan- 



tum Monte-Carlo algorithm 24 25 . The main challenge 



in this approach is to map the auxiliary systems back 
onto the real one, which can be accomplished by a stan- 
dard maximum entropy analytical continuation procedure 
26 . Details of the numerical procedure have been pro- 
vided in a recent publication 27 . Here, we compare the 
results of this Matsubara-voltage quantum Monte Carlo 
(MV-QMC) approach to data obtained with a scattering- 
states numerical renormalization group (SNRG) method 



29|30| and real-time quantum Monte Carlo (RT-QMC) 
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Model and Methods. — The simplest and most 
frequently used model for a quantum dot is the single- 
impurity Anderson model [3]. In this model, the dot is 
described by a single molecular orbital, which can ac- 
commodate up to two electrons. In the doubly occupied 
case the Coulomb repulsion leads to a charging energy 
of U . The dot orbital is connected via single-particle 
tunneling of amplitude ta to two continuous sets of non- 
interacting fermionic baths, which are called source and 
drain lead, and are indicated by an index a = ±1. The 
leads can have different chemical potentials, and the differ- 
ence e$ = — /i+i in chemical potentials represents the 
physical bias voltage. With these conventions, the Hamil- 
tonian reads 

H =^ eakaC'lk'^aka + ^ {ed+ (jB)d\da 



aka 



t . X 



aka ^ 



where cj^j.^. and are the usual fermionic creation op- 
erators of electrons with spin a = {±1} = {t?^}! in lead 
a = ±1 with momentum k or on the dot [h = c — = 1). 
The corresponding single-particle energies are eaka and 
e^, respectively, and we added the Zeeman energy aB = 
into the single-particle energy on the dot, which 
includes the effect of an external magnetic field H. We will 
consider the wide-band limit for the leads with symmet- 
ric coupling t- = t+ = t and the particle-hole symmetric 



point Cd = — J7/2. n denotes the phase space volume. As 
unit of energy we use the Anderson width F = 
where Mf is the conduction electron density of states at 
the Fermi energy. 

The method developed by Han and Heary employs a 
complexification of the physical bias voltage $ — > 'upm 
with iprn = 4:TTm/f3 (Matsubara voltages) 23 , which re- 
sults in an infinite set of auxiliary equilibrium statistical- 
mechanics systems one can efficiently solve by state-of-the- 



art Monte Carlo algorithms 24 25 . The mapping back to 



non-equilibrium quantities is done via analytical continu- 
ation i(pm — >■ $ ± i(5 from the Matsubara voltages to the 
real voltage. A central goal is of course the calculation 
of transport properties. Unfortunately, the non-locality of 
the current operator and the rather complicated analytical 
structure of the local Green's functions render such calcu- 
lations very difficult [28]. However, obtaining results for 
static local quantities, such as the double occupancy on 
the dot or the magnetization, is relatively straightforward 
27 . 

The full derivation of the formulas connecting those 
static quantities in the Matsubara voltage space with the 



real-voltage expectation value were presented in Ref. 27 
For local observables O, such as the double occupancy and 
magnetization of the dot, one finds the representation 



(0)(i^. 



{0)L 



■dip, (2) 



with go{^) the spectral function. The physical expecta- 



tion value is then given by 27 



(0)neq= {0)l 



VI '-^d^, (3) 



where Vj ■ ■ ■ denotes the principal value integral. 

Highly precise data for {0){iipm) can be obtained even 
for large values of ipm from the effective equilibrium sys- 
tems of the Matsubara voltage representation by use of the 
continuous-time quantum Monte Carlo (CT-QMC) tech- 
nique ^5| . In contrast to equilibrium QMC for the Ander- 
son impurity model, we expect a sign- respectively phase- 
problem here due to the presence of a complex quantity in 
the effective action. It turns out, however, that this phase 
problem is rather weak for small to intermediate bias, al- 
though it can become significant for large bias voltage and 
very large iyj^Q 

The representation ^ is formally similar to the stan- 
dard Lehmann-type representation of correlation func- 
tions. Due to the singular nature of the integral equa- 
tion ([2]), the numerical determination of the spectral func- 
tion Qoif) from QMC data is known to be an ill-posed 
problem. An adequate tool which helps to reduce un- 
controllable biases in the estimates of goiy^) is provided 
by the maximum entropy method (MaxEnt) 26 . It uses 



^For example, for iipm—60 we find for certain model parameters 
[U = ST, I3r = 20.0, B = O.OAT) average phase factors |{e'T)| 
0.073 at e$ = 2r, and |(e'T>| ^ 0.40 at e<S> = IT. 
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Fig. 1: Raw data used in the computation of the magnetization 
M and double occupancy D of a quantum dot {U = 8T, fiV = 
40, B = 0.02, e$ = O.Oir) in MV-QMC. The insets show the 
inferred offsets and spectral functions which yield the physical 
value. 



Bayesian inference by interpreting the spectral function as 
a probability density. A-priori and a-posteriori informa- 
tion about this density and quantities derived from it can 
be discriminated and discussed properly. As an illustra- 
tion of the procedure, we plot in Fig. [T] the raw data for 
the calculation of the magnetization M and double occu- 
pancy D in a model with U = 8T, l3T ^ 40, B ^ 0.02r, 
e$ = O.Oir. The offsets inferred from these QMC results 



are M. 



0.0643(1) and D 



= 0.1752(3) and 
the spectral functions obtained from the MaxEnt proce- 
dure are shown as insets. With these results, Eq. (|3| leads 
to M = 0.134(1) and D = 0.107(9). 

In order to analyze the numerical results of the 
Matsubara-voltage technique in the strongly correlated 
regime quantitatively, we compare them to SNRG [19(|29( 



30 and real RT-QMC 17 18 calculations 



The SNRG method extends the well-known equilibrium 
numerical renormalization group (NRG) [9] to describe 
steady-state non-equilibrium transport through interact- 
ing nano-devices. The methods starts from a formula- 
tion of the noninteracting U = Q problem in terms of 
exactly known scattering states, which are the solutions 
of the Lippmann-Schwinger equations 29 . Therefore, the 



boundary conditions of the open quantum system with 
particles entering and leaving the whole system are cor- 
rectly incorporated. Formulating the NRG in the basis 
of these scattering states allows for the application of the 
time-dependent NRG 31 32 . Starting from the initial 



U = Hamiltonian the interaction U is switched on and 
the time-evolution of the system in calculated. The re- 
sults for the long-time steady state limit can be accessed 
analytically and steady-state nonequilibrium expectation 
values are calculated for arbitrary interaction strength and 



bias voltage. 

The RT-QMC technique 17|18 is used to check the non- 
equilibrium double occupancies. This method computes 
steady-state expectation values by performing a quench, 
either in the interaction U or the voltage bias $, and is 
based on a stochastic sampling of weak-coupling diagrams 
within the Keldysh real-time Green's function approach. 
The technique employed here is the interaction quench, 
described in detail in Ref. [l8|. We start from the nonin- 
teracting system with applied bias voltage and switch on 
the interaction at time t — 0. The time-evolution of the 
double occupancy is then computed by randomly placing 
interaction vertices on the Keldysh contour, using a Monte 
Carlo technique. As a result of the quench, the double 
occupancy decreases and eventually approaches a time- 
independent value corresponding to the steady-state dou- 
ble occupancy of the interacting system 17 . If the times 



accessible in the RT-QMC calculation are long enough 
to see this convergence, the results are numerically ex- 
act. Difficulties arise in the small-voltage regime, where 
the transient dynamics becomes slow (the longest acces- 
sible time is limited by a dynamical sign problem in the 
Monte Carlo sampling). The double occupancy is easier 
to measure than the magnetic susceptibility, because in 
a half-filled dot with symmetric bias, only even pertur- 
bation orders contribute to the observable, which reduces 
the dynamical sign problem. 

Results. 

Double occupancy. The double occupancy D — (n-^nj^) 
is usually not discussed in the context of Kondo physics, 
and therefore we believe it is useful to start here with 
a study of the equilibrium behavior of this quantity in 
the absence of an external magnetic field, B = 0. We 
consider the particle-hole symmetric case where {nd) — 
{n^ + n^) = 1, and the double occupancy measures the 
charge fluctuation of the c^uantum dot, i.e. 

D = (ntn;) = \{{na if) = \{{na - {na)f). (4) 

This quantity is shown as a function of T in Fig. [2] for 
various values of the Coulomb repulsion U . The high 
temperature limit D{T F, [/) « (n^)(n^) = 1/4 is ap- 
proached by all curves (not explicitly shown) , and the low 
temperature saturation values roughly scale as 1/J7, as 
expected (see lower right inset of Fig. |2][a)). However, 
the first striking observation is that, in contrast to e.g. 
the magnetization, one does not see any explicit signature 
of the Kondo scale. Instead, for all values of U a mini- 
mum appears for temperatures on the scale of T ss F with 
weakly temperature dependent tails at lower T. The in- 
formation about Tk is contained in these tails. Plotting 
F{T) := 1 - D{T)/D(T = 0) as a function of T/Tk (main 
panel of Fig. [2]ja)) we find a nice scaling for T <T].^. The 
curves for all U fall on top of each other and follow a 
quadratic behavior, F(T — 0) ^ (T/Tk)"^, which is con- 
sistent with the Fermi liquid nature of the strong-coupling 
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Fig. 2; Equilibrium NRG results, (a) The scaling function 
F of the double occupancy D as a function of temperature in 
units of the Kondo temperature for various values of U . The 
dashed line shows the quadratic function O.OOla:^. The lower 
right inset shows the double occupancy D{T) as a function 
of temperature in units of F. The upper left inset shows the 
Kondo scale as extracted from a fit to the scaling function F ~ 
0.001{T/TKf as a function U. (b) The increase of the double 
occupancy relative to its value at the minimum, normalized to 
the total increase, Q^j.^gyi^Q'^y'"? ^ (Tmin is the temperature of 
the minimum). 

Kondo fixed point at T = 0. The Kondo scale, estimated 
from the onset of this scaling behavior via F{Tk) = 0.001, 
exhibits the expected exponential decrease as a function 
of U (see upper left inset). 

This unusual behavior must be interpreted in the fol- 
lowing way. For higher T, as U is increased, charge fluc- 
tuations are strongly suppressed by the Coulomb repul- 
sion and are frozen out as the system approaches its local- 
moment fixed point. This happens on a scale F, explaining 
why D{T) flattens on that temperature scale for all U. At 
very low temperatures, the systems tends to form a Kondo 
singlet due to spin fluctuations, which are accompanied 
by virtual charge fluctuations, thereby again enhancing 
D(T). (Note that the absolute value of this contribution 




Fig. 3: (color online) Double occupancy as a function of the 
bias voltage for a broad range of the bias voltage. The in- 
set shows the double occupancy for U = 5F, comparing to 
second-order perturbation theory (PT). At the present stage 
of development, significant mutual discrepancies are observed 
between all computational methods in the intermediate voltage 
range. 

to D[T) is of course again strongly suppressed with in- 
creasing U .) This increase produced by Kondo correla- 
tions should therefore occur at temperatures on the order 
of Tk. However, the enhancement of D{T) already set in 
for temperatures slightly below F ^ Tk which is observ- 
able in the lower right inset of Fig. [2](a) and Fig. [2]jb). 
This is due to the typical logarithmic tails [33] ubiquitous 
in Kondo systems as visible in Fig.[2](b). 

In Fig. [3] we provide non-equilibrium double occupancy 
data as a function of the bias voltage as obtained from 
the three approaches discussed above. The overall func- 
tional forms look quite similar to the equilibrium curves 
of Fig. [2] For large bias voltages all methods approach 
the non- interacting limit (n^n^) ^it^ = 1/4 and 

for zero bias voltage $ = they all reproduce the equi- 
librium value. (There are no RT-QMC data available for 
very small voltages, but it seems that the available data 
extrapolate to the SNRG and MV-QMC result.) At inter- 
mediate voltages, all methods produce a minimum. For 
the weakest interaction, U — 3F, all curves agree quanti- 
tatively. However, for larger U, the approaches differ both 
in the position and in the depth of this minimum. 

With increasing U, the MV-QMC approach exhibits the 
most pronounced minimum, whose position slightly moves 
to larger voltages. The RT-QMC results show a more 
shallow minimum, with a position roughly in agreement 
with MV-QMC. In contrast, the depth of the minimum in 
the SNRG results does not deepen significantly with U, 
while the position even appears to move slightly towards 
smaller voltages. 

In Ref. flSl it was shown that low order perturbation 
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theory results for the current are accurate up to U/T fa 4, 
with small deviations in the intermediate voltage regime 
for U/T = 6. At larger interactions, the perturbation 
theory result becomes incompatible with RT-QMC in the 
intermediate bias regime 1 < V/T < 5 [l8] . This is consis- 
tent with the findings of this study as shown in the inset 
of Fig. [3) where the RT-QMC results for U/T are close to 
the prediction from second-order perturbation theory. 

The origin of the discrepancies in the nonequilibrium 
double occupancy for intermediate and large U is unclear. 
If we assume that SNRG captures the correct position of 
the minimum in this would support the notion that 

a finite bias voltage can act similarly to an increase in tem- 
perature. The primary effect of both is to induce larger 
fluctuations, and therefore the behavior of the double oc- 
cupancy as a function of temperature and bias voltage is 
qualitatively the same. Consequently, the position of the 
minimum does not move significantly with bias voltage 
and temperature. 

On the other hand, the MV-QMC and RT-QMC resuhs 
suggest that voltage and temperature act in a fundamen- 
tally different way on the double occupancy. The shift of 
the minimum position to larger voltages with increasing U 
may be a genuine non-equilibrium effect, whose physical 
mechanism however remains to be clarified. 

Magnetic susceptibility. In Fig. [4j we present a com- 
parison of the SNRG and MV-QMC non-equilibrium data 
for the magnetic susceptibility as a function of the bias 
voltage for three magnetic fields B < Tk- The linear 
slope of the magnetization curves in a small, but finite 



magnetic field B, x 



is taken as an estimate 



for the susceptibility. Note that the computational effort 
for the MV-QMC calculations is kept constant through- 
out the considered parameter range, leading to a doubling 
of the relative statistical error as the magnetic field is de- 
creased by a factor of two. 

The magnetic susceptibility clearly exhibits a signature 
of the equilibrium Kondo scale, which is in contrast to 
the double occupancy. The susceptibility drops strongly 
when the voltage reaches the order of Tk indicating the 
destruction of the Kondo correlations due to the voltage 
induced fluctuations. As a function of temperature, this 
is a well-known fact in equilibrium [t], and can also be 
shown to hold for non-equilibrium systems, for example, 
using a simple scaling analysis of simulation data within 



the Matsubara-voltage formalism 27 



MV-QMC and SNRG agree reasonably well in the low- 
voltage regime. In both methods well-understood numeri- 
cal issues lead to a discrepancy with the equilibrium value 
of the magnetization at $ = 0. In MV-QMC, this is due to 
a systematic bias in the MaxEnt estimator and the growth 
in statistical error at low fields 27 . Within SNRG, the 



the main source of inaccuracies is the truncation of the ba- 
sis states due to the limited computer memory. Keeping 
the number of states fixed, this truncation becomes more 
severe in cases where many states become close in energy. 
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Fig. 4: (Color online) Magnetic susceptibility as a function 
of the bias voltage (in units of the Kondo temperature) ob- 
tained from MV-QMC and SNRG for [/ = SF, ^ = 40r-^ at 
different values of the magnetic field. Within MV-QMC the 
computational effort is chosen to be identical for each of the 
curves which renders B = O.OIF the least and B — 0.04r the 
most accurate. The equilibrium limit is shown as the red cross 
on the <3? = axis. The green dashed line shows a weighted 
least-square fit of the MV-QMC data to Eq. ([5}. We used 
the Haldane estimate of the equilibrium Kondo temperature 
Tk ~ ^F. The insets show close-ups of the low voltage region, 
as well as parameters of the fit to Eq. [5] 



which is the case when the voltage, the Kondo scale and 
the magnetic field are of the same order. A similar behav- 
ior of few-states SNRG can be observed for the equilibrium 
temperature dependence of m(T) (not shown) which can, 
however, be overcome computationally much more easily 
than in the non-equilibrium situation. 

While the sign problem in quantum Monte-Carlo pre- 
vents a direct computation of MV-QMC values at higher 



27 



voltages, this range is accessible with SNRG. In Rcf 
we noted that the Matsubara-voltage magnetization data 
can be described by the expression 



m(l>) 
B 



a 
B 



1 



((.2 



(5) 



where $ = <I>/(7rTK). The additional square-root in the 
denominator of Eq. ([5| leads to larger tails in the high- 
voltage regime of the susceptibility when compared to a 
fit with a pure Lorentz-like function. 



. This could 



be interpreted as a precursor of logarithmic tails expected 
in the Kondo regime 



33 



We used this function to fit the low-voltage MV-QMC 
data in Fig. |4] The first thing to note is that the low- 
voltage fit results in an excellent description of the SNRG 
data at high voltages, too. For B = O.OIF, where the 
MV-QMC data have the largest error, a fit with all three 
parameters in Eq. [5] became very bad. However, for larger 
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B with higher quaUty data, we observe that the param- 
eter c « 2. Constraining this parameter to c = 2 for 
B = O.Oir again results in an excellent fit here, too, in- 
cluding the high- voltage data from SNRG. Further analyz- 
ing the field dependence of the parameters a and 6, we find 
a surprisingly simple behavior, viz a « 0.75 • B/Tk and 
b K, 2{\ — B/Tk)- We therefore suggest as an approximate, 
but very accurate formula for the behavior of the local sus- 
ceptibility at low temperatures and fields T,B < Tk 



xiT,B,^) 



xiT,B,0) 



$2/4 



1 



(6) 



Summary. — In the present paper we have provided a 
comparison of three state-of-the-art computational meth- 
ods for two local observables on a quantum dot in a sta- 
tionary nonequilibrium state, namely the double occu- 
pancy and the magnetization. For the double occupancy, 
which only indirectly exhibits the Kondo temperature as a 
relevant energy scale, we have found substantial disagree- 
ment between all three methods, whose origin is unclear. 
However, a qualitative agreement between RT-QMC and 
MV-QMC is found, and for the interaction parameters 
considered, the RT-QMC gives results in close agreement 
with second-order perturbation theory. 

With regard to the magnetization, we compared MV- 
QMC and SNRG data within the Kondo regime. In con- 
trast to the double occupancy, we found good agreement 
between the two methods. This finding is remarkable, be- 
cause in these calculations, the physics is clearly controlled 
by the non-equilibrium Kondo effect. 
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